manuscript 



Dense packing crystal structures of physical tetrahedra 
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We present a method for discovering dense packings of general convex hard particles and apply it to 
study the dense packing behavior of a one-parameter family of particles with tetrahedral symmetry 
representing a deformation of the ideal mathematical tetrahedron into a less ideal, physical, tetra- 
hedron and all the way to the sphere. Thus, we also connect the two well studied problems of sphere 
packing and tetrahedron packing on a single axis. Our numerical results uncover a rich optimal- 
packing behavior, compared to that of other continuous families of particles previously studied. We 
present four structures as candidates for the optimal packing at different values of the parameter, 
providing an atlas of crystal structures which might be observed in systems of nano-particles with 
tetrahedral symmetry. 

PACS numbers: 61. 50. Ah, 45.70.-n, 02.70.-c 



Impenetrable (hard) mathematical bodies (e.g. 
spheres, spheroids, superballs, and polyhedra) have 
received much attention as models for the equilibrium 
behavior of systems of nano-particles and for the wealth 
of equilibrium and non-equilibrium structures they 
exhibit |lH6|. For the tetrahedron alone, quasicrystal 
structures, novel crystal structures, and glassy struc- 
tures have been reported in numerical simulations 
However, in all cases, the tetrahedral particles studied 
were mathematically ideal (polyhedral) tetrahedra. By 
contrast, in an experiment which found that regular 
tetrahedra have random packings that are denser than 
known for any other body, the tetrahedral macro- 
particles (dice) used had rounded edges and vertices 
PX)]. Tetrahedral nano-particles used as colloids are 
not only imperfectly-shaped tetrahedra, but are also 
soft, in the sense that the interactions beyond hard-core 
repulsion are significant [ll|, [l2|. In this paper, we 
attempt to characterize the packing behavior of physical, 
rather than mathematical, tetrahedra by studying 
a one-parameter family of particles with tetrahedral 
symmetry that interpolates between the mathematical 
tetrahedron on one end of the parameter's range and the 
sphere on the other end. We explore the effect of this 
parameter by constructing, using a de novo numerical 
search, candidate structures for the optimal packing of 
the different particles in the family. 

A particle interpolating between the sphere and the 
regular tetrahedron can be achieved by a variety of con- 
structions. The simplest construction, probably, is to 
place the centers of four unit spheres at the vertices of 
a regular tetrahedron with edges of length a and con- 
sider the volume at the intersection of all four spheres. 
We call the resulting figure a tetrahedral puff (Figure 
d]). For a special value of a, the tetrahedral puff is the 
Reuleaux tetrahedron (a three-dimensional version of the 
Reuleaux triangle, but not a solid of constant width) [l3|. 
A more convenient parameter than this edge length is 
the asphericity 7, which is the ratio between the radii 
of the particle's circumscribing and inscribed spheres. 
The value 7 = 1 obtained when the four spheres coin- 



FIG. 1: (Color online) tetrahedral puffs of varying asphericity. 
From left to right, the asphericities of the puffs shown are 
7 = 4/3,2,8/3. 



cide corresponds to a sphere. The value 7 = 8, which 
is the largest asphericity possible for a convex particle 
with tetrahedral symmetry and corresponds to a regular 
tetrahedron, is obtained in the limit that the four spheres 
intersect at a point. The Reuleaux tetrahedron is the puff 
with asphericity 7 = (3 + V24)/5 ^ 1.58. 

To efficiently search for candidate optimal packing 
structures we run a numerical search at increasing den- 
sities. To prevent overlaps between different particles 
in the structure we employ a recently developed over- 
lap resolution technique. In a packing of congruent 
particles, every particle can be obtained from a single 
primitive particle by a rotation and a translation: 
Ki = RiK + n = {Rix + ri|x G K}, where Ri, a ro- 
tation matrix, and ri, a translation vector, parameterize 
the configuration of the particle Ki . We call an exclusion 
projection the operation of, for any two particles Ki and 
K2 that overlap, identifying the new set of configuration 
parameters (R^, R2, r'l, ^^2) ^^^^ resolves the overlap while 
minimizing the distance in configuration-space from the 
original configuration, (Ri, R2, ri, r2). This is the projec- 
tion in configuration-space to the set of non-overlapping 
configurations. Below we give the method of implement- 
ing this projection. 

At high densities, however, we may have to resolve 
overlaps of a particle with multiple other particles, and 
the projection that applies to pairs cannot be used di- 
rectly. Therefore, we introduce multiple independent 
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copies (replicas) of the configuration parameters of each 
particle, so that after an exclusion projection, any pair 
of particles will have at least one overlap-free pair of 
replicas. Another projection ensures that in the final 
structure obtained from the search all replicas of a single 
particle agree on its configuration. This scheme, called 
divide and concur {D — C)^ has previously allowed us to 
study dense packings of a variety of polyhedral and high- 
dimensional spherical particles, and is here generalized 
to any convex particles. Apart from the exclusion pro- 
jection, described here, the rest of the application of the 
scheme to dense periodic packing discovery is described 
inRef. [14]. 

Given a pair of particles (or replicas) in a configuration 
(Ri, R2, ri, r2) which overlaps, the exclusion projection 
resolves the overlap by identifying a new set of config- 
uration parameters (R'^^, R2, r'^, r2) while minimizing the 
distance to the original configuration as defined by 

d2=^ ||r,-r:||2 + ||R,-R;||2, (i) 

i=l,2 

where the matrix norm is the Frobenius norm. The pro- 
jection algorithm for a general convex particle is most 
easily expressed in terms of the particle's support func- 
tion h{u) = maxxGK u • x. If the support function 
h{u) of K is known, then the support function of Ki 
is hi (u) = maxxGK u • R^x + u • = h{Rju) + u • [l5| . 

By the separating plane theorem, Ki and K2 do not 
overlap if and only if a vector u exists such that Ah{u) = 
hi{u) -\-h2{—u) < [HI. We can determine if such a vec- 
tor exists by numerically minimizing A/i(u)/| |u| |, which 
is bounded and attains a minimum over u. If the mini- 
mum value of A/i(u)/||u|| is positive, we must make the 
minimal change possible to the configuration parameters 
so that 

h[{u) + h^i-u) = for some u. (2) 

If we relax the condition that R^ is a rotation ma- 
trix, then we can reduce this constrained optimization 
problem to a simple unconstrained optimization problem 
in three vector variables. Namely, these vectors are u, 
vi = R']^u, and V2 = R2^u. Given these three vec- 
tors, the new configuration parameters which minimize 
([T]) and satisfy ([2]) are given by 

= ^1 - ^1^^^ ^^^^ 
^2 = r2 + ^^^^ (3b) 
R; = Ri + J\^\2(^f - ^^^^)' (3c) 

where Ah = h{vi) + /i(— V2) + (ri — r2) • u and the con- 
figuration distance is given by 

_2 A/1V2 + ||uRi - vilp + ||uR2 - V2IP ... 

^ = • (4) 



And so, we have reduced the problem, as promised, to an 
unconstrained minimization of @ over three vector vari- 
ables. Note that (|4]) is invariant under uniform positive 
rescaling of u, vi, and V2, and the resulting vanishing 
gradient direction must be taken into account when per- 
forming the minimization. 

The restriction on R^ to be a rotation matrix, i.e. the 
requirement on the rigidity of the particle, as in Ref. 
[1^ , is restored in the concurrence constraint of the D — 
C scheme. The projection of a general matrix into the 
subset of orthogonal matrices is as simple as taking its 
singular value decomposition and setting all the singular 
values to unity [16]. 

As a method for exploring dense configurations of 
general hard particles, we believe our projection-based 
method to be more direct and efficient when compared 
to event-driven MD simulations and stochastic MC meth- 
ods [2]. Partly, the D — C scheme draws its power from 
temporarily allowing non-physical configurations, with 
overlaps, non-concurring replicas, and non-rigid parti- 
cles, but then systematically acting to minimize these 
non-physicalities. In MC simulations, such temporary al- 
lowance has also been observed to be critical in efficiently 
exploring structures at high density [8]. 

For a selection of puffs of different asphericities, we use 
the D — C scheme to perform repeated de novo searches 
for periodic packings with p = 1, 2, 3, 4, 6, and 8 puffs per 
unit cell. The densities of the densest packings found for 
each value of 7 are given in Figure [3l Every packing den- 
sity that is reported here as putatively optimal for a puff 
of some asphericity has been reproduced by the numeri- 
cal search at least 3 times from random initial conditions 
and if the structure has two puffs in the primitive unit 
cell, it has been reproduced in searches with both p = 2 
and p = 4. 

As suggested by the results of our numerical searches, 
four different packing structures are optimal at differ- 
ent asphericities, separated by one continuous structural 
transition and two abrupt ones. Of most interest are the 
structures of the optimal packing for puffs of small as- 
phericity and large asphericity in the parameter ranges 
near the sphere and the tetrahedron respectively. The 
optimal packing structure for small asphericity, which 
we call the 5o-structure (for simple double lattice)^ is a 
tetragonal double lattice - that is, the union of two lat- 
tices (with one particle per unit cell) that are related to 
each other by an inversion about a point (Figure [2]) [17]. 
In the 5o-structure, the puffs are arranged into square 
layers so that from each layer, the puffs stick out on one 
side in parallel ridges running in one direction and on 
the other side in parallel ridges running in a perpendicu- 
lar direction. By stacking each consecutive layer with 
a 90° rotation, the ridges of one layer align with the 
ridges of the layer above it. In the limit 7^1, each 
layer approaches a square packing of spheres, and the 
5o-structure approaches the f.c.c. sphere packing struc- 
ture. 

For large asphericities, the optimal packing structure. 
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FIG. 2: (Color online) a unit cell of each of the four structures described, (from left to right) the So-, Pi-, <Si-, and Po -structures. 
For the So- and <Si-structure, a unit cell consisting of two primitive unit cells is shown. In all cases, the purple and pink puffs 
are related by inversion to the green and teal puffs respectively. The <So-structure is a body-centered tetragonal crystal where 
the body-centered puff is inverted in orientation from the corner puffs. The Pi -structure occurs when next-nearest square 
layers of the <So-structure come into contact, and its symmetry is broken by a re-orientation of different puffs of the same layer 
in different directions. By contrast, the <Si-structure arises by an abrupt transition at both ends of the parameter interval on 
which it is optimal. In the Po-structure the mirror planes of the green and teal puffs, which form a dimer, are not aligned, 
so that both can be in contact with the purple puff. In the limit 7^8, the green and teal puffs become aligned and form a 
bipyramidal dimer. The structure becomes the dimer double lattice structure reported in Refs. 0, 0] as the densest known 
packing of regular tetrahedra. 



which we call the Po-structure, is a dimer double lattice 
(i.e. each of the two inversion-related lattices is a lat- 
tice of dimers). The unit cell contains four puffs, two 
of which, which are in contact and form a shape simi- 
lar to a triangular bipyramid, are related by inversion to 
the other two (Figure [2]). We call each of the two in- 
version related pairs a dimer, in analogy with the dimer 
double lattice of Refs. [7, 9], which is the limit of the 
Po-structure as 7 ^ 3, and is the densest known pack- 
ing of regular tetrahedra. In this limit, each dimer ex- 
actly forms a triangular bipyramid. Note that away from 
7 = 3, the mirror planes of the two puffs constituting the 
dimer are not aligned with each other, making the dimer 
look twisted. This allows both puffs to form a contact 
with a nearby puff (see Figure [2|) , which in the limit of 
the bipyramidal dimer is facilitated simply by a contact 
along a common edge or vertex. 

The <So-structure appears to be the optimal packing 
structure from 7 = 1 to 7 ^ 1.63. On the other end of the 
asphericity scale, the Po-structure appears to be optimal 
from 7 ^2.19 to 7 = 3. However, in the intermediate 
range it appears that both of these structures are subop- 
timal and different structures take over. For 7 < 1.63, 
the density of the 5o-structure increases monotonically 
with asphericity. However, for the puff with 7 ^ 1.63, 
contacts between next-nearest layers of the <So-structure 
appear, and start to constrain the layer spacing. As- 
suming no change to the orientations of the puffs and to 
the construction of the layers, this constraint leads to a 
sharp drop in the density of the <So-structure. However, 
the structure found by the numerical searches shows a 



re-orientation of the puffs so that each layer is now com- 
posed of puffs of two different orientations (Figure [2]). 
This structure, dubbed the Pi-structure, still leads to a 
decline in the packing density, but a less dramatic one. 
Therefore, we have a local maximum in the packing den- 
sity at the transition from the 5o-structure to the Pi- 
structure. This transition seems to arise by a continuous 
deformation of the 5o-structure and is not abrupt. The 
re-orientation of the puffs suggests the beginning of a ten- 
dency towards the dimerization seen in the Po-structure. 
However, an intermediate structure is encountered be- 
tween the Pi-structure and the Po-structure, produc- 
ing another local maximum in the optimal density. This 
structure, to be called the <Si-structure, is a simple (non- 
dimer) double lattice without the tetragonal symmetries 
of the 5o-structure (Figure [2]), and is reminiscent of the 
simple double lattice structure reported in Ref. [7] . This 
structure appears to be separated from the others by an 
abrupt transition. Figure [3] plots out the densities of the 
different structures as obtained by the numerical searches 
and by analytic constructions. 

The proper context for the results obtained here for 
tetrahedral puffs is in comparison to two other one- 
parameter families of particles that include the sphere as 
a special case and whose dense packing structures have 
been investigated vigorously, namely spheroids [1] and 
superballs [2|] . The putative optimal packing of spheroids 
and superballs becomes monotonically denser the less 
sphere-like they become. By contrast, the optimal pack- 
ing density of puffs does not exhibit such monotonicity, 
although it is always higher than that of the sphere (con- 
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FIG. 3: Highest densities (filled markers) achieved for pack- 
ings of tetrahedral puffs of varying asphericity 7. The empty 
markers represent the density attained by a structure at an 
asphericity where it is suboptimal, as determined either by 
runs where the search got trapped away from the optimal 
structure or by runs where the number of particles per unit 
cell was incompatible with the optimal structure. The two 
abrupt structural transitions can be easily seen here where the 
density line for the <Si -structure crosses the lines for the other 
structures at critical asphericities at which the <Si -structure 
coexists with them. On the other hand, the continuous struc- 
tural transition between the So- and Pi -structures is associ- 
ated with a broken symmetry instead: the line corresponds 
to an analytic construction of the <So-structure, imposing its 
tetragonal symmetry. At 7 ~ 1.63, next-nearest layers come 
into contact with each other, leading to a sharp decline in 
density (dashed line). If we allow the tetragonal symmetry 
to be broken, but still allow only two particle orientations, 
the D — C search produces slightly higher packing densities 
(empty disks). The Pi -structure obtains higher densities by 
continuously breaking that symmetry as well. 



sistently with a conjecture by Ulam that the sphere is the 
worst-packing three-dimensional convex solid [18]). Un- 
like superballs, but like spheroids, the optimal packing 
of puffs is in all cases (besides the sphere) not a lattice 



packing, and the crystal unit cell includes at least two 
particles of different orientations. However, like super- 
balls, and unlike spheroids, the optimal packing struc- 
ture of puffs goes through an abrupt transition, where 
two dissimilar structures obtain an equal, optimal den- 
sity. 

A major difference of tetrahedral puffs in comparison 
to spheroids and superballs is the lack of inversion sym- 
metry. However, not only do the putative optimal pack- 
ings of puffs in all cases have such a symmetry (Figure 
[2]), it is in some cases the only symmetry of the pack- 
ing besides its lattice translations. Presumably, inver- 
sion symmetry plays an important role in forming close- 
packed structures of particles with tetrahedral symmetry 
and maybe even of other particles, a result already ob- 
served in the plane by Kuperberg and Kuperberg [17J. 

The tetrahedral puffs exhibit a much richer optimal 
packing behavior than either spheroids or superballs, and 
this richness is likely to be mirrored in the behavior of 
tetrahedral nano-particles. The variety of qualitatively 
different dense packing structures observed for mathe- 
matical tetrahedra is compounded when a physical shape 
parameter is added. A possible way to experimentally ac- 
cess the parameter investigated here, which describes a 
swollen tetrahedron, is by using colloidal particles that 
swell as a function of their temperature [l9|. Thus, a 
variety of structures and structural transitions could be 
explored. We have attempted here to provide an atlas of 
the possible crystal structures which might be observed 
in systems of particles with tetrahedral symmetry. Our 
candidates for optimal tetrahedral puff packings provide 
a starting point for the study of the phase behavior of sys- 
tems of particles with tetrahedral symmetry, both hard 
and soft, away from the limit of the mathematical tetra- 
hedron. The method presented here, applicable to any 
hard convex particle, could also be useful in characteriz- 
ing possible structures of many other particulate systems. 
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